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2 Abstract. We present results from global, three-dimensional magnetohydrodynamic 

3 (MHD) simulations of the solar wind/magneto sphere interaction. These MHD simulations 

4 are used to study ultra low frequency (ULF) pulsations in the Earth's magnetosphere driven 

5 by shear instabilities at the flanks of the magnetopause. We drive the simulations with 

6 idealized, constant solar wind input parameters, ensuring that any discrete ULF pulsations 

7 generated in the simulation magnetosphere are not due to fluctuations in the solar wind. The 

8 simulations presented in this study are driven by purely southward interplanetary magnetic 

9 field (IMF) conditions, changing only the solar wind driving velocity while holding all of the 

10 other solar wind input parameters constant. We find surface waves near the dawn and dusk 

11 flank magnetopause and show that these waves are generated by the Kelvin-Helmholtz (KH) 

12 instability. We also find that two KH modes are generated near the magnetopause boundary. 

13 One mode, the magnetopause KH mode, propagates tailward along the magnetopause 

H boundary. The other mode, the inner KH mode, propagates tailward along the inner edge of 

15 the boundary layer (IEBL). We find large vortical structures associated with the inner KH 

16 mode that are centered on the IEBL. The phase velocities, wavelengths, and frequencies of the 

17 two KH modes are computed. The KH waves are found to be fairly monochromatic with well 
is defined wavelengths. In addition, the inner and magnetopause KH modes are coupled and 

19 lead to a coupled oscillation of the low-latitude boundary layer. The boundary layer thickness, 

20 d, is computed and we find maximum wave growth for kd = 0.5-1.0, where k is the wave 

21 number, consistent with the linear theory of the KH instability. We comment briefly on the 

22 effectiveness of these KH waves in the energization and transport of radiation belt electrons. 



23 



1. Introduction 



24 One of the outstanding questions in the study of magneto spheric ultra low frequency 

25 (ULF) pulsations is the nature of their generation. Throughout this paper, when we refer to 

26 "ULF pulsations" we are referring to any broadba nd or quasi-monoch romatic pulsation in the 



27 range 0.5-15 mHz (Pc4-Pc5 bands, as defined by Uacobs et ail 111 96411 ). Several authors have 



28 shown that conditions in the sola r wind are well cor r elated with ULF pulsations observed in 



29 the magneto sphere. For example. \Mathie and Manm [|200jJ] show a strong correlation between 



30 solar wind speed and ULF pulsation power in the dayside magnetosphere, for L shells in 

31 the range L rs 4-7. The authors note that this high correlation is strong evidence that the 

32 Kelvin-Helmholtz (KH) i nstability at the magnetopause is the source of the pulsation energy. 



33 



Kepko and Spenca [1200311 conducted a study of a series of events in which ULF pulsations 



34 were observed in the dayside magnetosphere at a discrete set of frequencies. A spectral 

35 analysis of the solar wind density during the same time periods revealed significant wave 

36 power at the same set of discrete frequencies. This relationship suggests that variations in 

37 the solar wind dynamic pressure are responsible for driving ULF pulsations in the dayside 

38 magnetosphere. In addition, ULF variations in the Earth's convection electric field may 

39 respond dir ectly to variations in the o rientation and strength of the interplanetary magnetic 



40 field (IMF) Widley et all 



1997 



199811 



41 The suggested solar wind sources of magnetospheric ULF pulsations can be subdivided 

42 into three distinct driving mechanisms: pulsations observed near the dawn and dusk flank 



magnetopause driven by the strong velocity shear present there; pulsations in the dayside, 



44 driven by variations in the solar wind dynamic pressure; and pulsations driven by variations 

45 in the orientation and strength of the IMF. ULF pulsations generated by these different 

46 mechanisms are thought to occur primarily over different, but sometimes overlapping, loc al 



time sectors uakahashi and Anderson . 



19921 ; \Lessard et al. . 



1999; 



Ukhorskiy et all 



2005]. 



48 Thus, the global distribution of ULF wave power in the magnetosphere is an important 

49 diagnostic for understanding the generation mechanism(s). 



50 



The solar wind sources outlined above can be classified as external sources of ULF 



51 pulsations in the magnetosphere. In addition to these proposed external sources, a number of 

52 authors have suggested that processes internal to the magnetosphere may also be responsible 

53 for the generation of magneto spheric ULF pulsations. Wave particle interactions and local 

54 reconfigurations of the m agnetic field are b ut two examples of a number of proposed internal 



55 sources, see the review by 



Takahashi 



IU998I1 for more information. The focus of this paper will 

56 be on external driving of magneto spheric ULF pulsations and internally generated pulsations 

57 will not be discussed further. 

58 The spatial overlap of the distribution of ULF wave power for the different generation 

59 mechanisms complicates the study of the individual generation mechanisms. For example, it 

60 could be argued that a satellite measurement of a ULF pulsation in the dayside, near the dusk 

61 flank, was generated by either an impulsive variation in the solar wind density or driven by 

62 velocity shear, through the KH instability. Thus, a detailed knowledge of the upstream solar 

63 wind parameters is essential in determining the source of the ULF pulsation. This highlights 

64 one of the main difficulties in studying the three generation mechanisms proposed: there are 



65 very few events in which one of the three solar wind generating parameters is dominant over 



66 the other two. The solar wind is filled with complex structures and is quite dynamic. Typically 

67 all three of the suggested mechanisms are operating simultaneously. 
To circumvent these issues, we present results from a controlled experiment study of 



68 



69 ULF pulsations in the magneto sphere. We drive the Lyon-Fedder-Mobarry (LFM) global, 

70 three-dimensional, MHD simulation of the solar wind/magneto sphere interaction with 

71 idealized solar wind conditions. These idealized solar wind input parameters are chosen to 

72 mimic each of the three driving mechanisms outlined above. By holding all of the solar wind 

73 input parameters constant except one, we are able to study the effect of changing only that 

74 one parameter. The characteristics of the ULF pulsations generated by the particular driving 

75 mechanism under consideration can then be studied without the complications described 

76 above. The focus of this paper will be on ULF pulsations driven by the strong velocity shear 

77 near the dawn and dusk flank magnetopause. 
Magnetospheric ULF pulsations are also known to be im portant in the energization 



78 



79 and transport of radiation belt electrons. 



Rostoker et al. 



Ill998ll showed a strong correlation 



so between ou t er zon e electron flux and magnetospheric wave power in the ULF band. 



si \Baker et al.\ Ill998h similarly noted an association between ULF wave power and energetic 



82 electron enhancements in a comparison of two magnetic cloud events. For radiation belt 

83 electrons drifting in the equatorial plane, the most relevant field quantities for particle 

84 energization are the GSM z component of the magnetospheric magn etic field, B 7 , an d the 



as GSM azimuthal component of the magnetospheric electric field, jNorthrojx 1 1 96311 . Thus, 



se our efforts to characterize the ULF pulsations generated in the LFM simulations will be 



87 focused on pulsations in these two magnetospheric field components. Throughout this paper, 



88 we will comment on applications to radiation belt electron energization and transport, when 

89 appropriate. 

90 The remainder of this paper is structured as follows: In Section |2] we discuss the main 

91 theoretical and numerical work regarding the Kelvin-Helmholtz (KH) instability at the Earth's 

92 magnetopause. Section [3] provides a brief description of the global MHD simulations used 

93 in this study. In Section |4] we present the simulation results along with the spectral analysis 

94 techniques that are used to study the ULF waves in the simulation magneto sphere. Section [5] 

95 compares the simulation results with the theoretical KH results from Section [2] In Section [6] 

96 we provide a brief summary and concluding remarks. 

97 2. The Kelvin-Helmholtz Instability at the Magnetopause 

98 The Kelvin -Helmholtz instability o ccurs at t he int erface between two fluids in 



99 relative motion 



ChandrasekhaA 



1961]. 



Dungey 



195511 suggested that portions of the 



magnetopause boundary might be KH unstable. Observational evidence s uggesting a KH-typ e 



interaction at the magnetopause boundary soon followed. Surface waves mubry et al. 



02 \Levving and 



03 \Hones et al. 



1981 



Saunders et al. 



Burlaga, 19 791: \Fairfielai\l979c . \Sckopke et all 1 198111 and vortical structures 



1971 



198311 were observed propagating anti-sunward along the 



104 magnetopause boundary. 

105 Early th eoret i cal att e mpts to desc ribe t he KH inter a ction at the magnetopause boundary 



were done by \Sen\ n 1 96311 . weien ll 1 96411 . and lSouthwood IU968I1 . These linear MHD treatments 



all assumed the boundary interface between the magneto spheric and magnetosheath plasmas 



to be a tangential discontinuity (TD). A tangential discontinuity is a one dimensional layer 
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109 with velocities everywhere parallel to the planar interface. The total pressure and normal 

110 magnetic field are continuous across the interface. All three studies attempted to quantify 
m the effects of compressibility and found that for large relative flow velocities compressibility 
112 had a stabilizing effect. This is analog ous to hydrodynamic K H where it is well known that 



113 compressibility has a stabilizing effect nChandrasekhan\l96l\\ . This early work resulted in a 



necessary condition for the onset of the KH instability at the magnetopause 



115 is valid for incompressible plasmas separated by a tangential discontinuity nHasegawd 



boundary, which 



1975I1 : 



fk- v) 2 > 



(k-B 1 ) 2 + (k.B 2 ) 2 



(1) 



H rrii \ni n 2 J 

no B is the magnetic field, n is the number density, ji the permeability of free space, rrii the ion 

117 mass, k is the wave vector, and v the relative velocity between the two plasmas (v=Vi-v 2 ). 

us In Equation (OQ), the units are mks and the coordinate system is Cartesian with the boundary 

119 interface (e.g. the magnetopause) assumed to be planar. We define the boundary interface to 

120 be the YZ plane where the Y axis lies in the GSM equatorial plane, parallel to the boundary 

121 (positive tailward), the Z direction parallel to the GSM z direction and the X direction normal 

122 to the planar interface. Thus, the X and Y axes lie in the GSM equatorial plane, with the Y 

123 axis parallel to the boundary and the X axis normal to the boundary. In Equation (OQ), the 

124 subscripts 1 and 2 refer to the regions on either side of the planar interface, the YZ plane. 

125 We define X > to be region 1 and X < to be region 2. Along the dusk magnetopause, 

126 X > (region 1) corresponds to magnetosheath plasma and X < (region 2) corresponds 

127 to magneto spheric plasma. The wave vector k is restricted to the YZ plane (the boundary 
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128 interface). In what follows, we reserve capital XYZ for this boundary coordinate system and 

129 lowercase xyz for the standard GSM coordinate system used in our MHD simulations. Strictly 

130 speaking, Equation CD is only valid for incompressible plasmas separated by a tangential 

131 discontinuity; however, many features of the KH instability are well approximated in this limit 



133 



32 \Kivelson and Pii 



1984] 



The early theoretical KH treatments of \Sen\ II 196311 . \Fejen [1196411 . and \Southwood 



134 IU968I1 all assumed a tangential discontinuity at the boundary interface. However, satellite 

135 observations of magnetopause crossings revealed a thin ? viscous boundary layer at the 



36 magnetopause, dubbed the low-latitude boundary layer uiones et all 



1972 : 



Akasofu et al. 



37 \l973UEastman et a/.Ul976h . This boundary layer is roughly characterized by tailward flowing 



38 plasma on closed field lines. The existence of a thin boundary layer near the magnetopause 



139 suggested that modeling the magnetopause as a tangential discontinuity was inaccurate. In 

ho addition to this inaccuracy, a tangential discontinuity magnetopause cannot explain another 

hi key feature of observations: monochromatic surface waves. An incompressible KH model 

U2 that assumes a TD at the boundary interface predicts a growth rate (Equation ©, LHS-RHS) 

U3 that is a monotonically increasing function of the wave number, k. This implies a continuum 

144 of wavelengths will be excited and the smallest wavelength disturbances will grow the 

us fastest. This theoretical result contradicts magnetopause surface wave obser vations where 



46 monochromatic waves w ith well-defined wavelengths are typically seen [e.g. 



Takahashi et all 



1991 



Chen et al. . 



199311 



us The next level of sophistication in KH models came in the early 1980's where the effects 



U9 of compressibility and/or a boundary layer of finite thickness were included. The inclusion 



150 of either of these two effects complicated the calculations. Either the calculation of the 

151 characteristic equation remai ned analyt i cal but the roots, to (the comp lex frequency), had 



161 



163 



152 to be solved for numerically ttLee et al. 



1981 : 



Pu and Kivelson 



1983]. Or the linear MHD 



53 equat i ons were reduced to an eig envalue problem for u and integrated numerically Walker , 



54 ll98ll : lMara and PritchetA\19 82] 



The KH theory of Walken lll98lh included a boundary layer of finite thickness and 



156 assumed compressible plasmas. He showed that when the wavelength of the disturbance 

157 became comparable with the thickness of the boundary layer, the instability was quenched. 

158 This implied a fastest growing mode at a particular value of kd, where d is the boundary layer 

159 thickness. He studied the interaction for several geometric configurations and reported the 

160 fastest growing mode occurred for kd ~ 1. 



Results from a similar study (boundary layer/compressible plasmas) by \Miura and Pritchett 



62 IU982I1 f ound m a ximum wave growth for kd « 0.5-1.0 and were in good agreement with 



those of 



Walken Ill98lll . The reported values of kd at which maximum wave growth occurs 



164 should be interpreted qualitatively when applied to the real magnetopause. This is because 

165 the authors made various geometrical simplifications in their studies (B || v, B _L v) which 

166 are not always satisfied at the real magnetopause boundary. However, the main result from 

167 these two studies is clear: the KH instability will become quenched when the wavelength of 
lea the disturbance becomes comparable with the boundary layer thickness, i.e. when kd ~ 1. 

169 The value of kd at which the instability becomes quenched corresponds to the value of kd at 

170 which maximum wave growth will occur. Note that this result implies a particular wavelength 

171 for the fastest growing mode and thus, a particular frequency for the fastest growing mode (/ 
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172 = Vphase I A). The inclusion of a boundary layer of finite thickness is th us able to exp 



Walked lll98lh noted 



ain the 



173 observations of monochromatic waves with a well-defined wavelength. 

174 that the frequency of the fastest growing mode was in the Pc4-Pc5 range for typical values of 

175 k and d, inferred from observations. 

176 The inclusion of a bound ary layer of finit e thickness also allows for two KH modes to 



77 be generated at the boundary. \Lee et al.\ lll98lll included a boundary layer of finite thickness 



178 in their study of incompressible KH at the magnetopause. They reported that two KH modes 

179 were generated, one at the magnetopause boundary (the outer edge of the boundary layer) 

180 and one at the inner edge of the boundary layer (IEBL). They referred to these two modes 

181 as the magnetopause mode and the inner mode, respectively. They found the inner mode to 

182 be unstable most of the time whereas the excitation of the magnetopause mode depended 

183 critically on the orientation of the magnetic fie ld in the magnetosheath . It has been suggeste d 



184 that the vortical structures nHones et al. 



198 lh and the surface waves nCouzens et all 



198511 



as obse rved near the magnetopa use are associated with the KH instability at the IEBL. 



Pu and Kivelsom IU983I1 gave a comprehensive study of compressible KH at the 



187 magnetopause boundary. They assumed the boundary interface to be a tangential discontinuity 

188 and found two unstable KH modes, with different phase velocities and different wave vectors, 

189 k. They referred to these two modes as the fast and slow modes, where fast and slow 

190 refers to the different phase velocities. As with previous authors, they found the addition of 

191 compressibility to have a stabilizing effect. However, they found this effect to be small when 



192 compared with results in the incompressible lim it (Equation ©)• 



193 the apparent discrepancies in the early work of 



Senl Ill963h . 



heir treatment also resolved 



FeieA IU964I1 . and 



Southwood 



11 



195 



94 il968|l by recasting their results in terms of the s low and fast KH modes. 



A follow up paper 



96 [1 1 9 8 311 in the context of \Lee et at 



nKivelson and Pia 



1984] discussed the results of 



Pu and Kivelson 



Ill98lll . They noted that when the magnetopause and IEBL 



97 were separated by a large distance ( r elativ e to the amplitude of the disturbance) the fast 



98 and slow modes of \Pu and Kivelsom IU983I1 developed independently on the two interfaces. 



199 However, when the magnetopause and the IEBL were close together, the fast and slow 

200 modes coupled giving rise to two new modes, one mode propagating on the magnetopause 

201 (magnetopause mode) and the other propagating on the IEBL (inner mode). These two new 

202 modes had different phase velocities and different wavevectors, k. The phase velocity of the 

203 magnetopause mode was largely governed by the flow velocity in the magnetosheath while the 

204 inner mode phase velocity was governed by the flow velocity in the boundary layer. 

205 In what follows, we will compare results from global, three-dimensional MHD 

206 simulations of the solar wind/magneto sphere interaction with the theoretical results detailed 

207 above. We will demonstrate the existence of surface waves on the simulation magnetopause. 

208 These surface waves will be shown to be driven by strong velocity shear and not dynamic 

209 pressure variations in the solar wind. We will evaluate the condition for KH instability 

210 (Equation (OQ)) along the simulation magnetopause and show that it predicts the flow to be KH 

211 unstable at locations consistent with where the surface waves are seen in the simulation. We 



212 will use spectral analysis techniques to compute the frequency of these surface waves. We 

213 will also compute the wavelength of these surface waves directly from the simulation results. 
2u A simulation boundary layer thickness will be computed and the results will be shown to be 
215 consistent with kd = 0.5-1.0. We will also show that two KH modes are excited near the 
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216 simulation magnetopause boundary; one at the magnetopause and one at the inner edge of the 

217 boundary layer. We will present a scientific visualization of the simulation results that shows 

218 both of these KH modes propagating tailward along their respective boundaries. The scientific 

219 visualization will also reveal a coupled oscillation of the simulation boundary layer and large 

220 vortical structures associated with the inner KH mode. 

221 3. The LFM Global MHD Simulation 



222 The Lyon-Fedder-Mobarry (LFM) global, three-dimensional magneto spheric model 

223 solves the single fluid, ideal magnetohydrodynamic (MHD) equations to simulate the 

224 interaction between the coupled magnetosphere - ionosphere system an d the solar wind. T he 



225 details of the numerical methods used within the code are described in 



Lyon et al. 



[I2004h 



226 As an inner boundary condition, the magneto spheric portion of the code couples to a 2D 

227 ionospheric simulation which computes the cross polar cap potential, needed for the plasma 

228 flow boundary condition, based upon the field aligned currents at the inner spherical boundary 

229 and empirical models for the extreme ultraviolet and auroral conductances. The solar wind 

230 conditions, which form the outer boundary condition, can be taken from upstream satellite 

231 observations or can be created fro m scratch. Runs with r ealistic solar wi nd inputs have been 



232 used to study geomagnetic storms ^Goodrich et al. 



1998I1 and sub storms nLopez. et al. . 



19981. 



233 Idealized solar wind configurations have been particularly helpful in analyzing the physical 

234 processes involved in m agnetospheric phenomenon, such as the erosion of the magnetopause 



235 jWiltberger et all I2003I1 



236 



While the details of the numerical techniques used to solve the ideal MHD equations are 
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237 beyond the scope of this paper, they do have an impact on the simulations ability to resolve 

238 boundary layers. There are three key aspects of the numerical techniques used in the LFM 

239 that are important namely, the numerical order of the scheme, the use of nonlinear switches, 

240 and the size and shape of cells within the grid. The numerical order of a scheme can be 

241 thought of as the accuracy of the interpolation in terms of a Taylor series. A first order scheme 

242 introduces 'numerical' diffusion into the solution, while higher order schemes avoid diffusion 

243 at the cost of dispersion errors which introduce artificial extrema into the solution. Total 

244 variation diminishing (TVD) schemes are designed to balance the benefit s of high and fi rst 



245 order numerical schemes and are discussed in mor e detail in C hapter 21 of 



Hirsch 



[1988]. The 



246 LFM uses the Partial Donor Cell Method (PDM) nHairl I1987H as the nonlinear switch along 

247 with an eighth order interpolation scheme. In a simple test with linear advection problems, this 

248 approach allows for an increase by a factor of 400 in the Reynolds number when compared 

249 with a simple first order scheme. Since the numerical techniques used to solve the ideal MHD 

250 equations fall into the category of Finite Volume Methods, the cells used to discretize the 

251 computational domain are not required to be uniform or orthogonal. This allows us to place 

252 regions of high resolution in areas known a priori to be important, e.g. the magnetopause. 

253 In addition, these cells have aspect ratios designed to have more resolution in the directions 

254 transverse boundary than along it. In practice the numerical order and use of the PDM switch 

255 in the LFM are not changed, but we can adjust the grid resolution. In runs with the grid 

256 resolution changed by a factor of two in all directions we noticed roughly a 33% change in the 

257 thickness of the boundary layer. Simulations with another factor of two increase in resolution 

258 are not practical at this time. 
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259 To investigate the ULF pulsations generated by the strong velocity shear at the dawn and 

260 dusk magnetopause, we drive the LFM simulation with a range of idealized solar wind input 

261 parameters. The three LFM simulations used in this study differ only in the solar wind driving 

262 velocity. The remaining solar wind driving parameters are identical for the three simulation 

263 runs: B x = B y = nT, B z = -5 nT, n = 5 particles/cm 3 , v y = v z = km/s, and sound speed 

. 40 km/s. The three solar wind velocity inputs (corresponding to the three different simulations 

265 in this study) are v x = -400 km/s, v x = -600 km/s, and v x = -800 km/s. These idealized solar 

266 wind conditions are chosen to represent moderate driving of the magnetosphere system under 

267 3 different solar wind driving speeds. In order to allow the magnetosphere to take shape within 

268 the simulation domain, the IMF B z component begins with an interval of southward IMF, 

269 turns northward, and remains southward for the remainder of the simulation interval. The 

270 periods selected for analysis in this study are 4 hours long and occur two hours after the final 

271 southward turning of the IMF. The solar wind input parameters listed above are held constant 

272 during the selected 4 hours. Driving the simulations with constant solar wind parameters 

273 ensures that any discrete ULF pulsations in the simulation magnetosphere are not the result of 

274 perturbations in the solar wind. In particular, the solar wind dynamic pressure is held constant 

275 in these three simulation runs. Thus, any magnetopause surface waves that are generated 

276 cannot be the result of solar wind dynamic pressure fluctuations. From here on, we will refer 

277 to the three different simulation runs as the 400, 600, and 800 runs. 

278 The simulation results presented in this paper use a high resolution version of the 

279 magneto spheric grid. While the spacing between cells is not uniform in the region near the 

280 magnetopause, the typical cell size is approximately 0.125 Re (Earth radii). These simulations 
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are conducted with idealized solar wind conditions with no dip ole tilt in order to c oncentrate 



282 fully on the effects of velocity shear. As has been described by 



Korthetal 



ll2004h . the LFM 



283 does not produce significant region 2 field aligned currents or a ring current, which means that 

284 the fields in the inner magneto spheric portion of the simulation will be more dipolar than is 

285 seen observations. It also important to note that the LFM does not contain a plasmasphere and 

286 so the density profile in the inner magnetosphere will be different than the real magneto sphere. 

287 While these differences are important, they will not prevent us from examining the structure 

288 and evolution of magneto spheric ULF oscillations at the magnetopause flanks in a realistic 3D 

289 configuration. 



290 4. Simulation Results 

291 One of the advantages of this type of controlled parameter MHD study is the global, 

292 three-dimensional nature of the LFM MHD code. Analyzing the results from the three 

293 simulations provides a global picture of the distribution of ULF pulsations in the inner 

294 magnetosphere, under the three different solar wind driving speeds. We have developed 

295 a spectral analysis tool that provides a global map of where ULF pulsations occur in the 

296 simulation magnetosphere. We briefly describe this tool and the spectral analysis techniques 

297 used therein. 



298 



299 



300 



4.1. Spectral Analysis Techniques 

For the simulation field component of interest, say the simulation B z , we record a 4 hour 
time series at every spatial point in the simulation domain. At each spatial point, we compute 



301 



302 
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the one-sided, periodogram power spectral density estimate, P(f), of the zero-mean, 4 hour 
time series Xk, which we define as: 



P(fi) 



2dt 



X-j | 2 forj = 0,l,--- 



(2) 



N 



where 



N-l 



—27rijk 



N 



} for j = 0, 1, • ■ ■ , 



N-l 



k=0 



and 



J 



for J = 0, 1, • - ■ , — 



Ndt 



303 Here, dt is the sampling rate in seconds, fj are the discrete Fourier frequencies in Hz, N the 

304 number of points in the time series and X, the discrete Fourier transform (DFT) of the 

305 time series x/.. If the units of Xk are nT then the units of P(f) are (nT) 2 /Hz. For the three LFM 

306 simulations in this study, these parameters are dt = 30 seconds and N = 480. These sampling 

307 parameters determine the highest resolvable frequency, the Nyquist frequency, f^ y = 16.6667 

308 mHz and the frequency resolution, Af = 0.0694 mHz. 

309 The result of this computation gives P(f), the power spectral density estimate in the 

310 particular field component as a function of frequency, at every spatial point in the simulation 

311 domain. We can now build a global picture of ULF wave power in a given frequency band by 

312 computing, at each spatial point, the integrated power (IP) over a given frequency band of 

313 interest [/„, f b ], via Equation ©: 




fa 



(3) 
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3H which has units (nT) 2 in this example. Note t hat this quantity is different from the total po wer 



315 (TP) that is often used in ULF studies [e.g. 



Engebretson et al. 



1998 



Mathie and Mann . 



316 



20011: 



TP = Y J PU)) forall})e[f a ,f b } (4) 

3 

317 This quantity has units (nT) 2 /Hz in this example and should more accurately be called a total 

318 power spectral density. We favor Equation © over Equation © because Equation © does 

319 not explicitly account for the bandwidth, df. A better definition of TP would multiply the 

320 right hand side of Equation © by (f b — f a ) and thus, would have units (nT) 2 . Finally, we 

321 note that Parseval's theorem can be expressed in this terminology as the root integrated power 

322 (RIP) of P(f) equals root mean square (RMS) of the time series x^. 



\ 



k=0 



fNy 



P(f)df 



(5) 



(RMS = RIP) 



323 where f^ y = \j1dt is the Nyquist frequency. 

324 Computing power spectral densities from Equation © often results in noisy spectra when 

325 plotted versus frequency. Windowing the time series before computing the power spectral 

326 density estimate can smooth out this noisy behavior. When we need to examine the finer 

327 frequency details of our power spectra, we first window the time series with the discrete 

328 prolate spheroidal seq uences. This spectral estimation method is c ommonly referred to as the 



329 'multi-taper method' uhomson . 



19821 : 



Percival and Walden . 



199311 
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330 4.2. Spatial Distribution of ULF Wave Power 

331 Figure Q] shows the result of the prescribed technique for B z IP (top row) and E$ 

332 IP (bottom row) ULF wave power, integrated over the frequency band 0.5 to 15 mHz 

333 (Equation ©), for the 400, 600, and 800 km/s simulations (columns). Each panel is a GSM 

334 equatorial plane cut with 5 Re spaced ticks on the x and y axes (sun to the right). The black 

335 circle at the origin is the inner boundary of the simulation, located at r ~ 2.2 Re. The B z IP 

336 color scale ranges from to 75 nT 2 and the IP color scale ranges from to 5 (mV/m) 2 . 

337 The color scales in each row are the same to emphasize the increasing intensity of ULF wave 

338 power as the solar wind driving velocity is increased. The white contours in each of the 

339 panels in Figure [Hare B z =0 contour snapshots, which for these idealized solar wind driving 

340 conditions, is the approximate location of the magnetopause: The solar wind magnetic field 

341 is purely southward whereas the magneto spheric magnetic field is predominately northward. 

342 Thus, the B z =0 contour is good representation of the open/closed field line boundary. The 

343 bow shock is also resolved as the region of ULF wave power upstream of the B z =0 contour, 

344 particularly clear in the three B z IP panels (top row). 

345 Figure Q] shows substantial ULF wave power in the B z and field components near 

346 the dawn and dusk flank magnetopause. A close examination of the regions of intense ULF 

347 wave power shows that, in fact, there are three distinct ULF wave populations being driven 

348 in the simulations. The first distinction can be seen in Figure |2j which is taken from the 800 

349 km/s simulation. The leftmost panel in Figure [2] shows B z IP, integrated over the entire 

350 ULF band, 0.5-15 mHz (same panel as in Figured]). The middle panel in Figure |2] shows 
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351 B z IP, integrated over the frequency band 0.5-3 mHz. The far right panel shows B z IP, 

352 integrated over 3-15 mHz. The higher frequency population (3-15 mHz, Figure[2l far right 

353 panel) is confined to the magnetopause boundary whereas the lower frequency population 

354 (0.5-3 mHz, Figure [H middle panel) is interior the magnetosphere, away from the boundary. 

355 We will refer to the higher frequency (3-15 mHz) wave population generated near the 

356 magnetopause boundary as the Kelvin-Helmholtz (KH) population. In Figure [2l the two black 

357 ticks orthogonal to the magnetopause boundary mark the point on the magnetopause where 

358 the KH surface waves are first seen in the simulation. The lower frequency population (0.5-3 

359 mHz, Figure [2l middle panel) is generated by a process internal to the magnetosphere. This 

360 lower frequency ULF wave population along with its generation mechanism will be described 

361 in a follow up paper. In what follows, we will refer to this lower frequency population as the 

362 magneto spheric (MSP) population. 

363 The second distinction can be seen in the E$ IP panels (bottom row) in Figure [Hand is a 

364 distinction amongst the KH waves themselves. A close examination of the KH population near 

365 the dusk flank magnetopause in the 800 km/s, panel (bottom right) in Figure Q] reveals two 

366 distinct wave populations being driven near the dusk flank magnetopause (also true at dawn). 

367 In the panel, we see one region of intense ULF wave power aligned with the B z =0 contour and 

368 a second, spatially larger region of ULF wave power earthward of the magnetopause boundary. 

369 From here on, we will refer to the outer KH wave population, near the B z =0 contour, as the 

370 magnetopause KH mode and the more earthward KH wave population as the inner KH mode. 

371 We have verified that both the inner and magnetopause KH modes identified here in E^ IP 

372 are also identifiable in v r IP (not shown), which ensures that there are indeed two distinct KH 
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373 modes. 

374 To summarize, we have identified 3 distinct ULF wave populations driven in the 

375 simulations: the MSP population (Figured middle panel) and the two KH modes (Figured! 

376 bottom right), the magnetopause KH mode and the inner KH mode. This three-fold distinction 

377 is true for all three simulations in this study (the 400, 600 and 800 runs). 

378 4.3. Spectral Distribution of ULF Wave Power 

379 We now describe the spectral distribution of the ULF waves in frequency and wave 

380 number space. The spectral distribution of the ULF waves in frequency space is shown 

381 in Figure [3] Figure [3] shows radial profiles of B z (top row) and (bottom row) wave 

382 power spectral density (Equation ©) plotted along the dusk meridian for the 3 simulations 

383 (columns). The horizontal axis is distance along 1800 local time (LT), the vertical axis is 

384 frequency from to 14 mHz and wave power spectral density is plotted on the color scale. 

385 The vertical white lines represent the approximate location of the magnetopause. The ULF 

386 wave population excited near the magnetopause boundary is fairly monochromatic, with peak 

387 frequencies centered near 5, 8 and 10 mHz for the 400, 600 and 800 runs, respectively. The 

388 color scales are all different in Figure |3] so that the peak frequencies can be easily identified. 

389 The three panels in the bottom row of Figure [3] show both the magnetopause KH 

390 mode and the inner KH mode described in the previous section. The magnetopause KH mode 

391 is seen as the peaks in frequency centered on the white vertical lines (the approximate location 

392 of the magnetopause). The inner KH mode is seen as the peaks in frequency earthward of the 

393 white vertical lines. Note that the ULF wave power is more intense for the inner KH mode, 
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394 which can also be seen in the bottom row of Figure [IJ Also, the two KH modes have their 

395 peak power at the same frequencies, roughly 5, 8 and 10 mHz for the 400, 600 and 800 runs, 

396 respectively, which suggests that the two KH modes are coupled. 

397 The three panels in Figure [3] also show a limited radial penetration depth of the KH 

398 waves. The inner KH mode penetrates roughly 3 Re inwards of the magnetopause boundary, 

399 for each of the three solar wind driving velocities. This has implications for radiation belt 

400 transport an d energization where equatorially drifting el ectrons can be energized by these 



4oi ULF waves uiudson et at 



2000; 



Elkington et al 



2003I1 . However, as Figure [3] shows, this 



402 energization will only be effective within « 3 Re of the magnetopause boundary. This 

403 penetration depth is near the heart of the radiation belts (r ps 4-7 Re) only for the 800 km/s 

404 simulation, where the magnetosphere is highly compressed. However, the MSP population 

405 defined above (0.5-3 mHz, Figure [2l middle plot) is distributed rather uniformly along the 

406 entire dusk meridian, particularly clear in the B z panels in the top row of Figure [3j This 

407 population could effectively intera ct with radiation belt electrons through a drift resonant type 



408 interaction nElkington et all 1 199911 . 



409 An important quantity characterizing magnetospheric ULF pulsations is the azimuthal 

410 mode structure of the waves. Determining the azimuthal mode structure up to mode number 

411 m requires at least 2m simultaneous satellite measurements, distributed in azimuth. Thus, 

412 calculating the azimuthal mode structure from satellite measurements is especially difficult. 

413 Global MHD simulations are not limited by these criteria and are well- suited to study the 
4u azimuthal mode structure over a large range of m values. To calculate the a zimuthal mode 



415 structure, we follow the procedure outlined by jHokworth and Moz.eA 



1919] • This procedure 
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416 is essentially a Fourier transform in space followed by a Fourier transform in time. The spatial 

417 Fourier transform is done along circles of different radii. The result of the full procedure gives 

418 P(m, f), wave power spectral density as a function of azimuthal mode number and frequency, 

419 along different radii in the simulation domain. 

420 Figure |4] shows the E$ azimuthal mode structure for the 800 km/s simulation along three 

421 different radii in the simulation domain: 6.6, 8, and 10 Re- Here, and in Figure [3] above, the 

422 multi-taper spectral estimate described in Section |4~T1 has been used. In each of the color scale 

423 panels, the horizontal axis is azimuthal mode number from m = 0-30, the vertical axis is 

424 frequency from 0-15 mHz and E$ logarithmic power spectral density is on the color scale ( 

425 log(P(m, /)) ). The color scales are the same in each of the panels and range from -2 to 3. The 

426 bottom panels beneath each of the color scale panels show integrated wave power over three 

427 different frequency bands, 0.5-3 mHz (green), 3-15 mHz (red), and 0.5 -15 mHz (blue), to 

428 distinguish between the MSP and KH populations. The three panels in the figure show several 

429 interesting features. First, along the radius of 6.6 Re, we see the sub-3 mHz wave power, 

430 corresponding to the MSP population, and a hint of the KH populations near 10 mHz. As we 

431 move further out in radius to 8 and 10 R E , we begin to pick up the KH population near 10 

432 mHz. Second, the line plots underneath the three color plots show that the MSP and KH wave 

433 populations have their peak power at different azimuthal mode numbers. The MSP population 

434 (0.5-3 mHz, green) typically has its peak wave power near m ~ 8 and does not extend much 

435 beyond m ~ 15. On the other hand, the KH populations' (3-15 mHz, red) wave power is 

436 distributed over a much broader range of m values, say m rs 0-30, with its peak near m ~ 15. 

437 This feature is most evident in the 800 km/s, 10 Re panel (far right) where both populations 
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438 are being sampled. Similar features are seen in the azimuthal mode structure results from 

439 the 400 and 600 km/s simulations and in the B z azimuthal mode structure results from the 

440 400, 600 and 800 km/s simulations (not shown here): The MSP (0.5-3 mHz) population has 

441 its peak wave power near m ~ 8 and the KH populations' (3 -15 mHz) peak wave power is 

442 near m rs 15. 

443 It is interesting to note that the m number of the peak wave power for both the KH 

444 and MSP populations does not vary significantly as the solar wind driving speed is varied. 

445 This also has implications for radiation belt transport and energization where discrete peaks 

446 at a partic ular frequency and a particular mode numbe r will select the particles that will be 



447 energized ^Hudson et all 



2000; 



Elkington et ail 1 1 99911 . In particular, a given {m,f} pair 



448 will determine the drift frequency of the electrons that the KH waves could interact with, 

449 through the drift resonance condition: cu = m<jj d . Using this drift frequency, we can compute 

450 the relativistic first adiabatic invariant, M, for the given {m,f} pair. This value of M defines 

451 the particle population that could be energized by the KH waves. For this calculation, we 

452 use the dipole approximation for the L value and assume the electrons interact with the KH 

453 waves at the dusk meridian. Table [Hand Table |2] show the results of this calculation at two 

454 different points along the dusk meridian. Table \T\ shows the results for the most inward 

455 radial penetration of the inner KH mode. Table |2] shows the results of the calculation for the 

456 inner KH mode near the magnetopause. The values of the magnetic field, B, that are used 

457 in computing the relativistic correction factor are also shown. A 1 MeV electron drifting in 

458 the equatorial plane near geosynchronous orbit has an M value of roughly 1800 MeV/G. The 

459 values of M listed in Tableland Table[2]range from roughly 1/5 to 1/2 of this value. We thus 



460 



461 
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conclude that the KH waves could interact with equatorial plane electrons of a few hundred 
keV, near the L values listed in the table. 



462 5. Discussion 

463 5.1. Inner and Magnetopause KH Modes 

464 In order to fully characterize the distinction between the two KH modes alluded to 

465 above, we must first define a magnetopause boundary layer in the simulation. We define the 

466 simulation boundary layer (BL) as the continuous region of space that is 1 . earthward of the 

467 B z =0 contour (the magnetopause) and 2. where the local plasma flow is in the same sense 

468 as the local magnetosheath flow (tailward). Figure [5] is a GSM equatorial plane snapshot of 

469 the dusk flank magnetopause taken from the 400 km/s simulation (a scientific visualization 

470 of the simulation results can be downloaded here: ). Thi s scientific visualization was created 



471 with the CISM-DX visualization package for OpenDX nWiltberger et all 1200511 . The total 

472 electric field, |E|, is on the color scale, ranging from to 5 mV/m. We choose to plot |E| as 

473 opposed to Etj, because the two KH modes are most easily identified in |E|. The black vertical 

474 axis is the GSM positive y-axis, with ticks at 10 and 15 Re from bottom to top (sun to the 

475 right). The upper white contour is the B z =0 contour which is a very good approximation of 

476 the magnetopause in these idealized simulations. The lower white contour is a v x =0 contour. 

477 Near the dusk flank this contour tracks the approximate delineation between tailward flowing 

478 (boundary layer) plasma and non-tailward flowing (magnetospheric) plasma. The region 

479 between these two contours is approximately the simulation BL defined above. The black 
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480 contours are IP contours, outlining the inner KH mode and magnetopause KH mode 

481 populations (Figure Q3 bottom row). The two panels in the figure are identical except in the 

482 right panel we have replaced the IP contours with the local velocity field. 

483 The scientific visualization reveals the inner mode and magnetopause mode distinction 

484 described above. We see both the inner and magnetopause KH modes propagating along their 

485 respective boundaries. The inner mode is clearly seen as the tailward propagating blobs of |E| 

486 inside the larger black IP contour, just below the v x =0 contour. The magnetopause mode 

487 is less apparent. It propagates tailward inside the smaller black IP contour, just above 

488 the B z =0 contour. The coupled oscillation of the simulation boundary layer is striking. The 

489 st ructure of the s i mulat io n boundary layer is v ery similar to the diagram of Model B presented 



490 in \Sckopke et all IU981I1 . \Sckopke et al.\ IU981I1 proposed 3 models (A, B, and C) of the low 



491 latitude boundary layer to explain ISEE observations. Model A has both the magnetopause 

492 and the IEBL stable, Model B has both the magnetopause and the IEBL disturbed by surface 

493 waves and Model C has the magnetopause stable and the IEBL unstable. Our scientific 

494 visualization clearly shows both the magnetopause and the IEBL to be disturbed by surface 

495 waves and the BL configuration thus corresponds to Model B. A thickening of the simulation 

496 boundary layer through the KH region is also seen. The simulation boundary layer thickness 

497 near the right side of each panel in Figure [5]is roughly 0.5 Re and grows to roughly 1.3 Re 

498 near the left side of the panel. 



499 As discussed in the introduction to the \Safrankova et all Il2007h paper on variations in 



500 boundary layer thickness, there are many open questions reg arding the formation and structure 



of the low-latitude boundary layer. 



Song and Russeln II 199211 developed an explanation for the 
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502 formation of the LLBL during strongly northward IMF that relies on magnetic reconnection 



503 at high latitudes. 



Luhmann et al. 



IU984I1 presented a discussion for the formation of the 



504 bound ary layer on open field lines during southward IMF. Older work of 



Eastman and Hones 



505 Iil979ll indicated a role for viscous and diffusive mixing plasma from the magnetosheath onto 



506 closed field lines. In our results we are seeing antisunward flow on closed field lines during a 

507 prolonged interval of southward IMF, which we believe is reflection of the numerical viscosity. 

508 The right panel in Figure[5](and in the scientific visualization) shows the counterclockwise 

509 oriented vortices propagating tailward in the simulation BL. The orientation of the vortices 

510 is consistent w ith what is predicted by KH theory and w ith what has been observed near the 



511 magnetopause \Hones et alW\91u\Saunders et a/.LIl983ll . Note that the vortices are associated 



512 with the inner KH mode and are centered on the v x =0 contour, which is approximate 



1981 



Couzens et al 



y the 



198511 



513 IEBL. This fact has been alluded to many times [e.g. \Hones et al. 

514 Near the right side of each panel, where the KH waves are first seen in the simulation, a 

515 typical vortex size is roughly 1.7 R E in extent along the IEBL by roughly 1.0 Re in extent 

516 perpendicular to the IEBL. The vortices grow in size as they move downtail and can grow to 

517 be as large as roughly 5 R E by 3 R E near the left side of the panels. The ratio of the vortex 
sis dimensions in the equatorial plane remains constant at roughly 1.7 throughout the KH region. 

519 Kelvin-Helmholtz vortices are thought to be important for mass and momentum transport 

520 across the magnetopause, into the magneto sphere. This proposed mechanism is particularly 

521 important for northward IMF conditions when reconnection is less effective in plasma 



522 transport across the boundary jNykyri and OttoV l200ll ; \Hasegawa et all I2004I1 . It is certainly 



523 possible that the large vortical structures straddling the IEBL in the simulations could transport 
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524 plasma into the magneto sphere. However, in the present study, we make no attempt to quantify 

525 this possible transport mechanism. 

526 The coupled oscillatio n of the two KH m ode s is strong evidence that these are in fact the 



527 two KH modes described in 



Lee et al. 



IU981I1 and 



Kivelson and Pw II1984I1 . In order to confirm 



528 this, we must show the two modes have different phase velocities and different wave vectors, 

529 k. We can extract the phase velocity and wavelength characteristics of the two KH modes 

530 directly from the simulation results. By placing a line grid in the equatorial plane, through the 

531 two regions of E$ IP (Figure [6l left panel) and plotting along this line, we can calculate 

532 the wavelength for each of the modes. This corresponds to the wavelength in the Y direction 

533 in the boundary coordinate system defined in Section [2] The middle panel in Figure |6] shows 

534 this result for the inner KH mode in the 800 km/s simulation, from which we calculate a 

535 wavelength of Ay ~ 3.3 Re- The right panel in Figure [6]is essentially a time series of plots 

536 shown in the middle panel. Distance along the equatorial line grid is plotted on the horizontal 

537 axis, simulation time along the vertical axis and E$ is on the color scale. By measuring the 

538 slope of the linear features in the plot, we calculate a phase speed of fa 225 km/s. This panel 

539 also shows the coherent structure of the waves as they propagate downtail. Using Ay = 3.3 

540 Re and v p h aS e = 225 km/s, we calculate a wave frequency of / « v p h aS e/^Y ~ 11 mHz. This 

541 calculation of the wave frequency is in good agreement with the peak frequency observed in 

542 the far right panels in Figure [3] A similar calculation is done for the magnetopause mode in 

543 the 800 km/s simulation and for the inner and magnetopause modes in the 400 km/s and 600 

544 km/s simulations. The results are shown in Tableland Table SI Note that the Y direction is 

545 slightly different for the two KH modes. This is because the Y axis for each mode is chosen so 
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546 that it is parallel to the boundary for that mode; for the magnetopause KH mode, this boundary 

547 is the magnetopause and for the inner KH mode, this boundary is the IEBL. Thus, Ay should 

548 be interpreted as wavelength along each respective boundary. This slight difference can be 

549 seen in the left panel in Figure [6] as the two black lines (the respective Y axes) do not point in 

550 the same direction. 

551 The frequencies listed in Table [3] and Table |4] are in good agreement with the peak 

552 frequencies in Figure [3] This confirms that the two KH modes seen in the scientific 

553 visualization correspond to the two regions of E$ ULF wave power near the dusk flank 

554 magnetopause in the bottom row of Figure [IJ Moreover, when considering a particular 

555 simulation, the results in Table [3] and Table 0] show that the two KH modes have different 

556 phase velocities and different wavelengths but similar frequencies. For example, in the 800 

557 km/s simulation, we see that the phase velocity and wavelengths between the two KH modes 

558 differ by about 60 percent. However, the difference between the two frequencies is only about 

559 5 percent. A similar result holds for the 400 km/s and 600 km/s simulations. The coupled 

560 oscillation of the two KH modes is clear and we can positively identify t he two surface m odes 



561 in the simulation as the inner and magnetopause KH modes described 'm \Lee et al\ Ill98lll and 
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Kivelson and Pm IU984I1 



563 5.2. Boundary Layer Effects: Fastest Growing Mode 

564 The results from the previous section, along with the direct power spectral density 

565 computations (Figure [3]) show the two KH modes to be fairly monochromatic, with well- 

566 defined peak frequencies (Table [3] and Table ^Figure©. The monochromatic nature of the 
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567 waves is a direct result of the presence of a boundary layer, as discussed in Section [2] Recall 

568 that the KH instability is quenched when kd ~ 1 where k is the wave number and d is the 

569 boundary layer thickness. Thus, there is a particular wavelength (and frequency) for the fastest 

570 growing mode. We can therefore explain our discrete KH frequencies («5, 8, and 10 mHz 

571 for the 400, 600, and 800 km/s simulations, respectively) by showing that kd ~ 1 in our 

572 simulations. 

573 We begin by defining the wave appearance region (WAR) of the KH waves as the 

574 point along the magnetopause where the KH waves are first seen in the simulations. We 

575 determine these locations through a careful inspection of scientific visualizations from the 

576 three simulations. These points are located at 1624 LT along the magnetopause for the 800 

577 km/s simulation, at 1648 LT along the magnetopause for the 600 km/s simulation, and at 1708 

578 LT along the magnetopause for the 400 km/s simulation. For the 800 km/s simulation, this 

579 location is marked in the far right panel in Figure [2] with a black line perpendicular to the 

580 magnetopause. We use our simulation results near these points to calculate the simulation 

581 boundary layer thickness, as defined in Section I5T1 at the WAR. 

582 We compute the simulation boundary layer thickness near the WAR as follows: At the 

583 WAR, we extract the local velocity profile along a line perpendicular to the boundary (for 

584 example, Figure [2]). From this information, we compute the velocity locally parallel to both 

585 the magnetopause boundary and to the magnetosheath flow, in the equatorial plane. Figure [7] 

586 shows an example of this profile perpendicular to the magnetopause, at the WAR (1708 LT 

587 along the magnetopause), for a particular timestep in the 400 km/s simulation. The solid 

588 line is the parallel velocity plotted against distance orthogonal to the boundary. The vertical 
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589 dashed line indicates the point on the line orthogonal to the boundary where B z =0. This is the 

590 location of the magnetopause for this particular timestep. The vertical dotted line indicates 

591 the point on the line orthogonal to the boundary where the parallel velocity transitions from 

592 negative to positive values. This is the location of the IEBL for this particular timestep. The 

593 distance between these two vertical lines is the simulation boundary layer thickness, d, at the 

594 WAR. In Tableland Tabled we show the results of this computation for the boundary layer 

595 thickness, d, at the dusk WAR, for the three simulations in this study. We note that near the 

596 WAR, the simulation boundary layer thickness fluctuates throughout the 4 hour interval. The 

597 values of d listed in Tableland Table [6] are the average values for the 4 hours of simulation 

598 time and are typical values for the thickness depth. It is also important to note that the LFM 

599 grid resolution near these points is sufficient to resolve this boundary layer thickness. There 

600 are typically 3-4 grid cells within the simulation boundary layer. 

601 In order to evaluate kd, we must also calculate k. In Equation (QQ), the wave vector k 

602 is restricted to the YZ plane. We can approximate k from our computed values of Ay listed 

603 in Tableland Table |4] (i.e. k ~ ky). This is a reasonable approximation, as can be seen 

604 in Figure [8] where is plotted on the colorscale from -6 to 6 mV/m in the YZ plane for 

605 the inner KH mode in the 800 km/s simulation. As described in Section |2l the Y axis lies 

606 in the equatorial plane and is parallel to the boundary, in this case the IEBL. The Z axis is 

607 parallel to the GSM z axis. The origin of the coordinate system in Figure [8] is located on the 
60s magnetopause at the WAR (1624 LT). The axes ticks are spaced at 1 Re and a black line 

609 that makes a 20° angle with the equatorial plane is also shown. Note that the KH waves are 

610 generated near the equatorial plane, which can be inferred through a careful inspection of 
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en Equation Clearly, the KH waves propagate not only in the positive Y direction (tailward) 

612 but also in the Z direction. This indicates a small kz component to k, in addition to ky = 2n 

613 /Ay and that the approximation k pa k Y is valid. The results of the kd calculation are shown 
6u in Tableland Tabled under the assumption k pa ky. Our values of kd are consistent with 

615 kd ~ 1. This explains why we see KH waves of a particular wavelength (or frequency) in the 

616 simulation results. The presence of the boundary layer of finite thickness quenches the KH 

617 instability and thus we have maximum growth for a particular k and a particular frequency, /. 

618 The monochromatic KH waves seen in the simulations are manifestations of this process. 

619 In Section I5T1 we compared the frequencies computed directly from the simulation 

620 results (Table |3] and Table H]) with the peak frequencies from the power spectral density 

621 computations. Similarly, we can compare the peaks in azimuthal mode number for the KH 

622 modes (Figure with the wave numbers computed directly from the simulation results 

623 (Tableland Table©. In order to do so, we must transform the k Y values computed in the 

624 boundary coordinate systems into the GSM coordinate system where the azimuthal mode 

625 structure calculations were done. Thus, we must simply decompose ky into k r and k^ = m/r. 

626 The results of this decomposition are shown in Table [7] and Table [8] for the two KH modes. 

627 We see that the azimuthal mode number, m, lies between 12 and 19 for both KH modes and 

628 all three solar wind driving speeds. These values of m are in good agreement with the peaks 

629 in power spectral density seen in Figure |H where we found m pa 15 for both KH modes and 

630 all three solar wind driving velocities. For a particular simulation, the values of m listed in 

631 Tableland Table[8]show a slight difference in m between the two KH modes. This difference 

632 cannot be resolved from the power spectral density computations shown in Figure |4] due to the 
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633 narrow azimuthal width separating the two KH modes. 

634 5.3. Criteria For KH Instability 

635 From the simulation results, we can directly evaluate the condition for KH instability 

636 (Equation (OQ)) to see where it predicts the flow to be KH unstable. As Equation (OQ) is only 

637 valid for a tangential discontinuity, we make no attempt to evaluate it in the simulation 

638 boundary layer. For this calculation, we assume that there is no boundary layer and use the 

639 field values on either side of the boundary layer, outside of the boundary layer. For example, 

640 for region 2 (the magneto sphere) fields, we use field values that are earthward of the IEBL. 

641 Similarly, for region 1 (the magnetosheath) we use fields that are away from the magnetopause 

642 and in the magnetosheath proper. Equation (OQ) cannot predict whether the inner KH mode or 

643 the magnetopause KH mode or both are excited. It can only predict whether the field values in 

644 the magnetosheath proper and the magnetosphere proper are such that the KH instability will 

645 or will not occur. There is only one KH mode in the incompressible, tangential discontinuity 

646 KH theory that is used to derive Equation ©. 

647 All of the field quantities in Equation (OQ) are specified by the simulation results. For the 

648 wave vector k, we use the ky values listed in Table [5] for the inner KH mode. We choose 

649 the inn er mode ky valu es as the inner mode is predicted to be the more unstable of the two 



650 modes nLee et q/.Ul981ll . We evaluate this condition along the equatorial plane magnetopause, 

651 from subsolar past the dusk flank, and we assume that k is parallel to v. This is a reasonable 

652 assumption given that the calculation is done in the equatorial plane and that the fastest 

653 growing mode will occur for this orientation of k and v. The results of this calculation are 
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654 shown in Figure [9] for a particular timestep in the 800 km/s simulation. The horizontal axis 

655 is LT along the magnetopause and the vertical axis is the left-hand side (LHS) minus the 
ese right-hand side (RHS) in Equation (QQ). The horizontal dashed line corresponds to marginal 
657 stability. The vertical dashed line marks the point along the magnetopause where the KH 

ess surface waves are first seen in the simulation, i.e. the WAR, as defined above. For the 800 km/s 

659 simulation, this point is located at 1624 LT along the magnetopause. The trace of LHS-RHS 

660 shows that the condition for KH instability is first satisfied somewhere near 1400 LT along 

661 the magnetopause. We now address the question of why the KH waves are not seen in the 

662 simulation until points near 1624 LT on the magnetopause; Figure [9] suggests they should first 

663 appear somewhere near 1400 LT. 

664 We begin by noting that a positive value of LHS-RHS in Equation (QQ) is the square of the 

665 linear growth rate of the KH waves. Thus, Figure [9] shows the square of the linear growth rate 

666 of the KH waves as a function of distance along the equatorial magnetopause. In Figure [9l 

667 we see that near 1400 LT, where the condition for KH instability is first met, the square of the 
ees growth rate is rs 0.0055, so that the growth rate is ps 0.0742 in this region. Thus, the e-folding 

669 time in this region is ps 2ir I 0.0742 = 85 seconds. We can now calculate the growth length 

670 in the region between 1400-1624 LT from this e-folding time and an estimation of the phase 

671 speed near 1400 LT. A plot of the magnetosheath speed parallel to the magnetopause (not 

672 shown here) shows the value of the magnetosheath flow speed to be ps 260 km/s near 1400 LT 



673 Thus, the value of the KH phase velocity in this region is ps 260/2 km/s =130 km/s Walker , 



674 Il98lh . These two calculations imply that the growth length in the region between 1400-1624 



675 LT is « 130 km/s * 85 s = 1.7 Re- Thus, the waves will travel along the magnetopause 
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676 a distance of roughly 1.7 Re from 1400 LT before they can grow to a sufficient size to be 

677 resolved in the simulation. Finally, we note that the magnetopause arc length between 1400 

678 and 1624 LT is roughly 5.4 R E . This partially explains why the KH waves are not seen in the 

679 simulation until points near 1624 LT along the magnetopause. The waves do not grow to a 

680 resolvable size until they travel roughly 1.7 Re along the magnetopause. We now calculate 

681 an improved estimate of the growth length based on a more applicable KH theory in order to 

682 explain the disparity between the growth length of 1.7 Re predicted by Equation £0 and the 

683 value of 5.4 R E . 

684 The disparity between where the waves are seen in the simulations and the growth 

ess length calculation done above is probably due to the unrealistic assumptions used in deriving 

686 Equation (OQ). Equation (OQ) is valid for incompressible plasmas separated by a tangential 

687 discontinuity. The LFM simulation solves the compressible MHD equation s and the resolves 



ess a realistic magnetopause boundary layer. The KH theory of 



Walker] [| 198 ill , which solves 



689 the compressible MHD equations in the presence of a boundary laye r, is a more ac curate 



ego description of the KH instability at the magnetopause. In particular, vNalken lll98lh finds 

691 maximum (normalized) wave growth rates for ^DjV in the range 0.1-0.3, where 7 is the 

692 growth rate, D is half the boundary layer thickness, and V Q is half the relative velocity between 

693 the two plasmas. Near 1400 LT, where Equation (OQ) first predicts the flow to be KH unstable, 

694 the value of jV is roughly 0.9, using the 7 value near 1400 LT (0.0742), and the simulation 

695 values near 1400 LT for D (3121/2 km) and V Q (130 km/s). Thus, Equatio n (Q]) predicts a 



696 normalized growth rate that is much larger than what is reported in 



Walked 11198 lh . Assuming 



697 a normalized growth rate of -yD/V = 0.25 and using the LFM simulation results near 1400 LT 
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698 for D and V Q , we calculate a growth length of 6.2 Re- This calculation of the growth length 

699 is in better agreement with the distance between where Equation <0Q) first predicts the flow to 

700 be KH unstable and where the waves are first seen in the simulation, a distance of roughly 5.4 

701 Re- Similar results hold for the 400 km/s and 600 km/s simulations (not shown here). 

702 At this point, it should be clear that the surface waves seen near the dawn and dusk flanks 

703 in the three simulations are indeed Kelvin-Helmholtz waves. The simulation surface wave 



704 



characteristics are consistent with the theoretical and observational KH surface wave results. 



705 The simulated waves have the proper frequencies, wavelengths, phase velocities, propagation 

706 directions and they have the large vortical structures associated with them. Furthermore, we 

707 see maximum wave growth for values of kd consistent with theoretical predictions. We also 
70s find that the theoretical results predict the magnetopause boundary to be KH unstable and the 

709 theoretical growth rate of the waves is consistent with where the waves are first seen in the 

710 simulations. Again, we emphasize the fact that the solar wind dynamic pressure is constant 

711 in our simulations. Thus, the surface waves cannot be attributed to fluctuations in the solar 

712 wind dynamic pressure, a cla im that is often us ed to discount observational evidence of KH 



713 generated surface waves \e.g. \Song et a/.LIl988ll 



7H As an aside, we note that the KH instability has been invoked to explain surface 
715 waves and vortical structures se en in global MHD simulations driven by real solar wind 



716 conditions. 



Slinker et al. 



Il2003h compared LFM simulation results with Geotail observations 



717 of magnetopause crossings. The LFM simulation reproduced the surface waves observed by 

718 Geotail a nd the authors noted that th e likely source of the oscillations was the KH instability. 



719 Similarly, 



Collado-Vega et al. 



Il2007h simulated 9 hours of a high speed solar wind stream that 
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720 was seen at LI between 29 March to 5 April 2002, using the LFM simulation. The authors 

721 reported large vortical structures near the magnetopause boundary and attributed these vortices 

722 to the KH instability. In both of these studies, the authors suggested the the KH instability was 

723 responsible for the surface waves and vortical structures but offered no conclusive evidence 

724 that the KH instability was indeed the source. 

725 Finally, we note that all of the KH theory discussed in this paper is linear MHD wave 

726 theory. Thus, once the KH waves have developed into their nonlinear stage, the linear wave 

727 theory is no longer applicable. The formation of the large vo rtical s t ructures in the si mulation 



731 



1984 



Wu 



1986D. Thus, 



728 is strong evidence that we have reached the nonlinear stage nMiurd 

729 applying the linear theory at points along the magnetopause boundary where the waves have 

730 reached their nonlinear stage is invalid. 



6. Summary and Conclusions 



732 In this paper, global, three-dimensional MHD simulations of the solar wind/magneto- 

733 sphere interaction were used to study ULF pulsations in the inner magneto sphere. The 

734 MHD simulations were driven with idealized, constant solar wind input parameters. These 

735 parameters were chosen to study the effect of changing only the solar wind driving velocity, 

736 while holding the other solar wind input parameters constant. Driving the simulations with 

737 constant solar wind parameters ensured that any discrete ULF pulsations in the simulation 

738 magnetosphere were not driven by fluctuations in the solar wind. The simulation results 

739 revealed ULF surface waves near the dawn and dusk flank magnetopause. These surface 

740 waves were shown to be driven by the Kelvin-Helmholtz instability and not dynamic pressure 
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741 fluctuations in the solar wind. 

742 A closer examination of the surface waves revealed that two KH modes were seen near 

743 the dawn and dusk flank magnetopause. These two KH m odes were identified as the inner 



744 KH m ode and the magnetopause KH mode, as described m \Lee et al.\ ll 198 lh : \Kivelson and Pu 



745 II1984I1 . The magnetopause KH mode was found to propagate tailward along the magnetopause 



746 boundary whereas the inner KH mode was found to propagate tailward along the inner edge 

747 of the boundary layer (IEBL). These two KH modes were found to have different phase 

748 velocities and different wavelengths but oscillated at the same frequency. We presented a 

749 scientific visualization that showed the coupled oscillation of the two KH modes and a coupled 

750 oscillation of the low-latitude boundary layer. The scientific visualization also revealed large 

751 vortical structures associated with the inner KH mode. These vortical structures were centered 

752 on the IEBL and propagated tailward along the IEBL, growing in size as they moved downtail. 

753 Both KH modes were found to occur for kd = 0.5-1.0 where k is the wave number and d is 

754 the boundary layer thickness. This fact was used to explain the monochromatic nature of the 

755 KH waves. The frequency of the KH waves was found to depend on the solar wind driving 

756 velocity, with larger driving velocities generating KH waves with higher frequencies. The 

757 azimuthal mode number, m, of the KH waves was found to be between 15-20 and did not 

758 change significantly with solar wind driving speed. The relativistic first adiabatic invariant, 

759 M, was computed from the m and / values of these KH waves. We found that the KH waves 

760 could effectively interact with equatorial plane radiation belt electrons of a few hundred keV, 

761 near the dusk meridian. 



Figure 1. 
Figure 2. 
Figure 3. 
Figure 4. 
Figure 5. 
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913 Figure Captions 




Figure 1. Global distribution of ULF integrated power (IP, 0.5-15 mHz) in the GSM equa- 
torial plane for the simulation B z (top row) and (bottom row) field components. The three 
columns correspond to the three MHD simulations used in this study (v sw = 400 km/s, 600 
km/s and 800 km/s, respectively). The white contours are B z =0 contours, the approximate 
location of the magnetopause. The KH surface waves are manifest as the regions of intense 
IP near the dawn and dusk flank magnetopause (sun to the right, 5 Re spaced ticks). The E^ 
IP panels in the bottom row show the two distinct KH populations, the inner KH mode and 
magnetopause KH mode. The color scales in each row are set to the same value to emphasize 
the increasing intensity of ULF wave power as the solar wind driving speed is increased. 
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Figure 2. B z IP in the GSM equatorial plane from the 800 km/s simulation, integrated over 
three different frequency bands to highlight the Kelvin-Helmholtz (KH) and magnetospheric 
(MSP) ULF wave populations. The left panel is B z IP integrated over 0.5-15 mHz (same 
panel as in Figured)). The middle panel is B z IP integrated over 0.5-3 mHz to highlight the 
MSP population. The right panel is B z IP integrated over 3-15 mHz to highlight the KH 
population. In each panel, the two black lines perpendicular to the magnetopause mark the 
point along the magnetopause where the KH waves are first seen in the simulation. 
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Figure 3. Radial profiles of B z (top row) and (bottom row) power spectral density along 
the dusk meridian for the three simulations in this study (columns). Distance along 18LT is 
on the horizontal axis, frequency is on the vertical axis and power spectral density is on the 
color scale. The vertical white lines represent the approximate location of the magnetopause. 
The KH waves excited near the magnetopause boundary are fairly monochromatic with peak 
frequencies near 5, 8, and 10 mHz for the 400, 600 and 800 km/s simulations, respectively. 
The three panels in the bottom row show both the magnetopause KH mode (peaks in fre- 
quency near the magnetopause) and the inner KH mode (peaks in frequency earthward of the 
magnetopause). Note the limited radial penetration depth of the inner KH mode (bottom row) 
and the uniform distribution of the MSP population (0.5-3 mHz) across a substantial portion 
of the dusk meridian (top row). 
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Figure 4. azimuthal mode structure in the 800 km/s simulation, along three different radii 
in the simulation domain: 6.6, 8, and 10 Re- In the color scale panels, the horizontal axis 
is azimuthal mode number, m, the vertical axis is frequency and logarithmic power spectral 
density is on the color scale. The three panels beneath each of the color scale panels show 
integrated wave power (IP) over three different frequency bands, 0.5-3 mHz (green), 3-15 
mHz (red), and 0.5 -15 mHz (blue), to distinguish between the MSP (green) and KH (red) 
populations. Note that the MSP population has its peak wave power near m ~ 8 whereas the 
KH population has its peak wave power near m ~ 15 (far right panel). The same is true for the 
400 km/s and 600 km/s simulations (not shown here). The color scale is the same in all three 



panels and ranges from -2 to 3. 
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Figure 5. Scientific visualization snapshot of the dusk flank magnetopause from the 400 km/s 
simulation. The black vertical axis is the GSM positive y axis with ticks at 10 and 15 Re 
(sun to the right). In both panels, |E| is on the color scale from to 5 mV/m. The upper 
white contour is a B z =0 contour, the approximate location of the magnetopause. The lower 
white contour is a v x =0 contour, the approximate location of the IEBL. The region in between 
these two contours is the simulation boundary layer. In the left panel, the two black E$ IP 
contours are shown to outline the inner and magnetopause KH modes. The inner KH mode 
propagates tailward inside the larger black IP contour, near the IEBL. The magnetopause 
KH mode propagates tailward inside the smaller black E^ IP contour, near the magnetopause. 
In the right panel, the E^ IP contours are replaced with the local velocity field. Note that the 
counterclockwise oriented vortices are associated with the inner KH mode and centered on the 
IEBL. These vortices grow in size as they propagate tailward from roughly 1.7 Re by 1.0 Re 
near the right side of the panel to roughly 5 R E by 3 R E near the left side of the panel. Also 
note that the boundary layer thickens through the KH region, from roughly 0.5 Re near the 
right side of the panel to roughly 1.3 Re near the left side of the panel. 
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Figure 6. Example of how the wavelength and phase velocities are computed for the two KH 
modes in the simulations. The left panel is a dusk flank zoom in of the IP panel in FigureQ] 
for the 800 km/s simulation. The two black lines in this panel define to the Y directions for 
the two KH modes. The middle panel shows plotted along the black line for the inner KH 
mode, from which we calculate a wavelength of Ay « 3.3 Re- The right panel is a time series 
of plots shown in the middle panel. The horizontal axis is distance along the inner mode black 
line (left panel), the vertical axis is simulation time and E^ is on the colorscale. By measuring 
the slope of the linear features in this panel, we calculate a phase velocity of v p hase ~ 225 km/s. 
Note the coherent structure of the waves as they propagate downtail. 
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Figure 7. Parallel velocity profile near the KH wave appearance region (WAR) for a particular 
timestep in the 400 km/s simulation. The velocity parallel to both the magnetopause bound- 
ary and to the magnetosheath flow (i>||) is plotted along a line perpendicular to the boundary 
(horizontal axis). The vertical dotted line is the location of the IEBL while the vertical dashed 
line is the location of the magnetopause. The region between these two lines is the simulation 
boundary layer, as defined in the text. We see a boundary layer thickness, d, of roughly 0.65 
Re- 
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Figure 8. The extent of the KH waves out of the equatorial plane, for the inner KH mode 
in the 800 km/s simulation. The origin of the coordinate system is located at the KH wave 
appearance region (WAR), 1624 LT along the magnetopause. The Y axis lies along the IEBL 
in the equatorial plane with the positive direction tailward. The Z axis is parallel to the GSM 
z axis. is on the color scale from -6 to 6 mV/m. The axes ticks are spaced at 1 Re- Note 
that the KH waves are generated near the equatorial plane and propagate in both the Y and Z 
directions. 
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Figure 9. Evaluation of the condition for KH instability (Equation (OQ), k 1 1 v) along the equato- 
rial plane dusk magnetopause for the 800 km/s simulation. The solid trace shows left-hand side 
(LHS) minus right-hand side (RHS) from Equation ©. The horizontal dashed line corresponds 
to marginal stability. The vertical dashed line marks the point along the dusk magnetopause 
where the surface waves are first seen in the simulation (1624 LT, the WAR). Note that the con- 
dition for KH instability is first satisfied somewhere near 1400 LT. In Section 1531 we explain 
why the KH waves are not seen in the simulation until points near 1624 LT. 
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9H Tables 





{m,f} 


L'min 


(nT) 


M 
(MeV/G) 


400 Run 


{15,5} 


10.5 


31 


377 


600 Run 


{15,8} 


8.6 


52 


469 


800 Run 


{15,10} 


7.5 


79 


513 



Table 1. Relativistic first adiabatic invariant M values computed for the given {m,f} pair, for 
the three simulations in this study. These M values determine the electron populations that the 
KH waves could interact with. We also show the L and B values along the dusk meridian where 
we assume the interaction occurs. These values are for the most inward radial penetration of 
the inner KH mode. 





{m,f} 


Lmax 


B(L max ) 
(nT) 


M 
(MeV/G) 


400 Run 


{15,5} 


13.5 


15 


601 


600 Run 


{15,8} 


11.6 


16 


738 


800 Run 


{15,10} 


10.5 


24 


835 



Table 2. Same as Tabled] except the M values are computed near the magnetopause 
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400 km/s 


600 km/s 


800 km/s 


Vphase (km/s) 


140 


160 


225 


Ay (R E ) 


4.2 


3.3 


3.3 


/(mHz) 


5.2 


7.6 


10.7 



Table 3. Inner KH mode equatorial plane phase velocities and wavelengths, computed directly 
from the simulation results, and the resulting wave frequencies. 





400 km/s 


600 km/s 


800 km/s 


Vphase (km/s) 


180 


225 


375 


Ay (R E ) 


5.2 


4.3 


5.2 


/(mHz) 


5.4 


8.2 


11.3 



Table 4. Same as Table[3]for the magnetopause KH mode. 





400 km/s 


600 km/s 


800 km/s 


d(R E ) 


0.53 


0.48 


0.47 


k Y (I/Re) 


1.50 


1.90 


1.90 


kd 


0.80 


0.91 


0.90 



Table 5. Simulation boundary layer thickness, d, the Y component of the wave vector k in 
the boundary coordinate system, and the product kd, for the inner KH mode in the three sim- 
ulations (under the assumption k « k Y \ see Figure [8]). Note the values of kd in the range 
0.5-1.0. 
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400 km/s 


600 km/s 


800 km/s 


d{R E ) 


0.53 


0.48 


0.47 


k Y (I/Re) 


1.21 


1.46 


1.21 


kd 


0.64 


0.70 


0.57 



Table 6. Same as Table[5]for the magnetopause KH mode. 





400 km/s 


600 km/s 


800 km/s 


k r (URe) 


0.68 


0.86 


1.04 


kff, (I/Re) 1 rn 


1.33/18 


1.69/19 


1.6/16 



Table 7. The equatorial plane components of the wave vector k in the GSM coordinate system 
for the inner KH mode in the three simulations. Note that the azimuthal mode number, m, lies 
between 12 and 19 for all three solar wind driving speeds, in good agreement with the m peaks 
in Figure HI 





400 km/s 


600 km/s 


800 km/s 


k r (11 Re) 


0.66 


0.79 


0.70 


kfj, (I/Re) 1 rn 


1.02/16 


1.22/17 


0.99 / 12 



Table 8. Same as Table[7]for the magnetopause KH mode. 



